Mon. Not. R. Astron. Soc. 000, 000-000 (1997) Printed 1 February 2008 



Stability of Power-Law Disks I — The Fredholm Integral 
Equation 

N. W. Evans and J. C. A. Read 

Theoretical Physics, Department of Physics, 1 Keble Rd, Oxford, OXl 3NP 
1 February 2008 

ABSTRACT 

The power-law disks are a family of infinitesinially thin, axisymmetric stellar disks 
of infinite extent. The rotation curve can be rising, falling or flat. The self-consistent 
power-law disks are scale-free, so that all physical quantities vary as a power of radius. 
They possess simple equilibrium distribution functions depending on the two classical 
integrals, energy and angular momentum. While maintaining the scale-free equilibrium 
force law, the power-law disks can be transformed into cut-out disks by preventing 
stars close to the origin (and sometimes also at large radii) from participating in any 
disturbance. 

This paper derives the homogeneous Fredholm integral equation for the in-plane 
normal modes in the self-consistent and the cut-out power-law disks. This is done by lin- 
earising the coUisionlcss Boltzmann equation to find the response density corresponding 
to any imposed density and potential. The normal modes - that is, the self-consistent 
modes of oscillation - are found by requiring the imposed density to equal the response 
density. In practice, this scheme is implemented in Fourier space, by decomposing both 
imposed and response densities in logarithmic spirals. The Fredholm integral equation 
then relates the transform of the imposed density to the transform of the response den- 
sity. Numerical strategies to solve the integral equation and to isolate the growth rates 
and the pattern speeds of the normal modes are discussed. 
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1 INTRODUCTION 



This paper begins an investigation into the large-scale, global, linear modes of a family of horizontally hot, vertically 
cold, idealised disk galaxies with rising, falling or flat rotation curves. Here, we collect all the mathematical and 
numerical details needed for our study. A companion paper completes the investigation by presenting results on the 
global spiral modes, together with the astrophysical implications. The focus is exclusively on the in-plane instabilities, 
such as the modes causing bar-like or lop-sided distortions. Of course, our hot, razor-thin models are also unstable 
to bending modes (Merritt & Sellwood 1994). These are ignored because they almost certainly can be eliminated by 
giving the disks a modest thickness. 

In our study, we have the good fortune to be able to follow a magisterial earlier analysis carried out by 



Zang (197(:). This remarkable Ph. D. thesis was supervised by Alar Toomre and also benefited from a number 
of ingenious suggestions provided by Agris Kalnajs. What Zang did in the mid-seventies was to carry out the first 
complete, global, stability analysis for a family of differentially rotating, stellar dynamical disks. The object of his 
study was an infinitesimally thin disk in which the circular velocity is completely flat. The special case in which 



the stars move on circular orbits is known as the Mestel (1963) disk. Their hot stellar dynamical counterparts are 



usually called the Toomre-Zang disks. The global spiral modes in this model were described extensively by Zang. 

The power-law disks are infinitesimally thin disks in which the circular velocity varies as a power of cylindrical 
polar radius, viz., Vdvc oc R~^^'^. They are amongst th e simplest ste llar dynamical disks known - the distributions 



of stellar velocities that build the models are given by Evans (1994) . Their rotation curves can be rising or falling 



The self-similarity of these disks simplifies the analysis considerably. It enables much of the global normal mode 
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calculation to be performed exactly. The model examined in Zang's (1976) thesis is the particular case with a 
completely flat rotation curve (/3 = 0). The perfectly self-similar disks are somewhat special. To supplement our 
analysis, we have also concentrated on a modified version of the pure power-law disks, in which the central regions 
of the disks are immobilised or cut-out. This was achieved by reducing the fraction of active stars - those able to 
respond to any disturbance - from unity at moderate radii to zero at the centre. In order to aid comparison of 
our results with numerical simulations, especially those of Earn (1993), we have also tapered the fraction of active 
stars to zero at large radii. The immobile components still contribute to the potential experienced by the stars. The 
material in Zang's (1976) thesis wa s never com pletely published, although some of these and various subsequent 
calculations are reported in Toomre (1977; 198l| ). So, we urge the reader to remember that our contribution consists 
only in generalising the analysis of Zang from the disk with a completely flat rotation curve to the entire family of 
power-law disks. This is not trivial, as the limit of a flat rotation curve is a singular one. Nonetheless, our job has 
been very considerably eased by being able to lean on the sturdy support of Zang and his two implicit co-workers. 
We use the notation Z followed by an equation number as a convenient shorthand for reference to Zang's (1976) 
thesis. 

The aim of this paper is to derive the homogeneous, singular, Fredholm integral equation for the normal modes 
in the self-consistent and the cut-out power-law disks. First, the equilibrium properties of the models are introduced 
in Section 2. The coUisionless Boltzmann equation is linearised to derive the Fredholm integral equation in Section 3. 
Section 4 summarises the strategies for its computational solution. This paper describes only methods. The following 
paper in this issue of Monthly Notices presents the numerical results and astrophysical consequences. 



2 THE EQUILIBRIUM DISKS 

This section begins with a brief introduction to the power-law disks. Section 2.2 introduces new variables - the 
home radius and the eccentric velocity - which can be used to characterise the stellar orbits according to their size 
and shape. The frequencies and the periods of orbits are derived in the next section. The equilibrium distribution 
functions of stellar velocities of the self-consistent power-law disks are presented in Section 2.4. The concluding 
section explains how to "cut out" the centre of the disk and to taper the disk at large radii. Both these require 
changes to the distribution function. 



2.1 The Potential-Density Pair 

The equilibrium density of the power-law disks is 



Scq - So ( 



(1) 



The self-consistent potential in the plane of the disk is (3chmitz fc Ebert 1987; Lemos et al. 1991; Evans 1994) 
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where the reference velocity vp is defined as 
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Clearly (0) fails for /? = 0. In this case, we have (Mestel 1963) 



(2) 



(3) 



(4) 



Whenever an expression fails at /3 = 0, the corresponding result for the Toomre-Zang disk should be used. It can 
often be derived by taking the limit (3 ^ i) with rHopital's theorem. The circular velocity is 



Up. 



(5) 



Thus, the reference velocity vp is the circular velocity at the reference radius Rq. Disks with /3 > have falling 
rotation curves, whereas disks with /3 < have rising rotation curves. Disk galaxies typically have flatfish rotation 
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curves (e.g. Rubin et al. 1978 ; [Mihalas fc Binney 1987 ) - sometimes slowly rising, sometimes slowly falling at large 
radii ( |Casertano fc van Gorkom 1991 ). The total mass of the disk is infinite for all /3 in the range — 1 < /3 < 1. The 
mass enclosed within a radius R is 



M{R) 



1-/3 



R_ 

Rq 



1-/3 



(6) 



2.2 The Home Radius and the Eccentric Velocity 

The Lagrangian for stars orbiting in a general power-law disk is 



2 / D \ /5 



There are two isolating integrals of motion, namely the energy E and the angular momentum L^. In terms of the 
radial velocity u — R and the tangential velocity v ~ R9, the integrals are 



Following Zang ( 1976| ), let us define the home radius Rn to be the radius at which the tangential velocity is equal 
to the circular velocity. By conservation of angular momentum, we have 



Again following Zang (1976), the eccentric velocity U is defined as the maximum radial speed reached during an 
orbit. U is thus positive by definition. It is easy to show that the eccentric velocity is attained at the home radius 
(^ and that its value is given by 

— (I-O(f)"' 

A similar derivation using the potential of the Mestel disk gives the result already well-known to Zang (Z2.29) 

= 2E-v^(l + 2ln-^Y (11) 



VqRq 



We shall find it convenient to work in dimensionless co-ordinates. We define the following dimensionless integrals of 
motion: 



v} [RoJ ' i?o' [RoJ ' ^ vpRo- 



Likewise, dimensionless radial and tangential velocities, radius and time coordinates can be defined as: 

The scaled energy (8), home radius (^) and eccentric velocity ( p^ are given by 

E^^{u'+v^)--^, Rh^LI^, [72 = 2^-1 + |. (14) 

We will also need expressions for the scaled radial and tangential velocities: 

2 /~a .\ - 1 



|(«--l), (15, 
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In terms of the scaled coordinates, the equations of motion derived from the Lagrangian become 

d^R _ I 1 

Imagine solving the equations of motions with starting positions and velocities. If two stars have the same u and 
V, then their orbits have the same shape, although the size of the orbit is proportional to i?H- This leads us to an 
intuitive understanding of the eccentric velocity and the home radius. Stellar orbits with the same dimensionless 
eccentric velocity have the same shape. The dimensionless home radius fixes the overall size of the orbit. 

The star's orbit in the equilibrium disk is limited by its isolating integrals of motion. In terms of the dimensionless 
co-ordinates, the turning-points of the orbit occur when 

= U^ + 1- R-'^ + ^ {r-'^ - l) = 0. (17) 

For given U and /3, there are two solutions, i?min corresponding to pericentre, and i?max corresponding to apocentre. 
These are easy to find numerically by the Newton-Raphson technique. The star thus moves within an annulus of 
the disk. For certain values of the eccentric velocity, the orbit closes, and the star then traverses a one-dimensional 
manifold within the annulus; in general, the orbit does not close, and the star eventually passes arbitrarily close to 
every region of the annulus. 



di 



1 

R2- 



(16) 



2.3 Periods and Frequencies 

Let us start by defining a useful auxiliary integral 



Jn{U) - 2 / r, (18) 

jRr..AU) ^„ (^72 + 1 - i?-2 + I - 1^ 



where -Rmin and iJmax are the solutions of (|l_3)- For the [3 — Q case, this becomes (Z2.39) 

Jn{U) = 2 I ^ -. (19) 
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To evaluate J'n numerically, we remove the singularities at either end of the integrand by transforming to a variable 
9. We define R = m + asinO, where m is the midpoint of the radial motion, m — ^(i?min + -Rmax) and a is its 
amplitude, a — i(i?inin — ^max)- The integration can then be carried out using the midpoint method (Press ei 



al. 1989, chap. 4). 



The radial period T is the time taken for the star to travel between two successive pericentres: i.e., to move out 
and in again. The corresponding radial frequency is defined by k = 2tt/T. Using the symmetry of the orbit, we have 

T = 2 r""'^^ dt = 2 dR^2_R^(R^^^ f-- d_R ^^^^ 

Jr=R,„.„ jR„.,n R '"P \RoJ Jr^„,„ u 



Using the expression for the radial velocity given in (15), along with the definition (18), we find that the radial 
period and frequency are 

We shall find it useful to define the dimensionless radial frequency k: 

vp \Rq J Jq 

The angular period Q is defined to be the angle through which a star moves during the time taken to complete one 
radial oscillation: 

e = = 2 . 2 = 2 1™ = J.. ,23) 

Jt=o yfl.„.„ R Jr^,„ uR uR 
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Figure 1. The radial frequency ft and angular frequency f2 plotted against eccentric velocity U, for /3 = -0.9, -0.6, -0.3, 0, -1-0.3, +0.6, 
+0.9. In the left-hand plot, the dotted lines on the vertical axis mark the epicyclic limit of ^2 — fi. As U — > oo, the frequencies k and 
remain positive but become vanishingly small. 



The angular frequency is defined to be the mean angular speed of the star: = 8/T. We shall commonly use the 
dimensionless angular frequency (l, where 

fi 

The dependence of k and fl on U and /3 is shown in Fig. ^. In the limit [/ — s- 0, the dimensionless epicyclic and 
circular frequencies are 



Ko = lim k= yj2 - P, no = lim = 1. (25) 

In the limit C/ — > oo, these frequencies tend to zero from above. As expected, the angular frequency is that of a star in 
a circular orbit at i?H, namely flo = Vch-c/Rn- At low eccentric velocity, there is superimposed on this circular motion 
a small radial oscillation with frequency kq = V2 — P I'circ/^H- This is of course just the epicyclic approximation for 



near-circular orbits (Binney & Tremaine 1987) 



2.4 The Structure of Phase Space 

The distribution function of the stars can depend on positions and velocities only through the isolating integrals 
of motion, according to Jeans' ( 



1919D theorem. Even quite simple disks can possess rather complicated distribution 



functions (Evans & CoUett 1993). The power-law disks are attractive candidates for a stability analysis because they 
possess a rich family of simple self-consistent distribution functions fs found by Evans ( |1994| ) . These are built from 
powers of energy and angular momentum, and normalised so that the integral of the distribution function over all 
velocities recovers the surface density. They are: 



where 
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and 



7 > -1; 



7>-/3-l. 



(27) 



(28) 



Note that these formulae differ by a factor of two from those given in Evans' (1994) paper. This is because Evans' re- 
sults assume a bi-directional disk, where the stars rotate in both senses. In this paper, we consider only uni-directional 
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disks. Analogous distribution functions for the Toomre-Zang disk have been known for a long time (Bisnovatyi- 
Kogai 1976; Zang 1976; Toomre 1977; Binney & Trcmaine 1987, chap. 4) 



/,(£;,L,) = (7L2exp 
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where C 
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In these formulae, 7 is a constant anisotropy parameter. Its physical meaning will shortly become apparent on 
examining the dynamical quantities derived from the distribution function. 
The mean streaming velocity, or the stellar rotation curve, is 
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(30) 



(31) 



(32) 



The radial velocity dispersion a^i is the root-mean-squared radial velocity of the stars. We shall also find it convenient 
to define a dimensionless radial velocity dispersion du- These are: 



1 + 7 -f 2/3 V i? 
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(33) 



If all the stars move on circular orbits, the radial velocity dispersion vanishes. Such a disk is said to be "cold" , by 
analogy with the motion of molecules in gases. As the eccentricity of the orbits increases, the stars acquire random 
motions and thus the "temperature" of the disk increases. High values of the anisotropy parameter 7 correspond to 
low velocity dispersions, i.e., cold disks. This is obvious on examining the form of the distribution function: fs cc LJ. 
When 7 is large, more of the stars have high angular momentum and lie on near-circular orbits. The isotropic model 
is given by 7 = 0. The squared tangential velocity dispersion is the difference between the second tangential 
velocity moment (w^) and the square of the mean streaming velocity (v)'^. The second moment (w^) is derived from 
the distribution function as: 



1 + 7 + 2/3 



Ro 
R 
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(34) 



As 7 c», crfj and (w^) — > Wcirc, so the disk is rotationally supported. These distribution functions have the 
property that, at any spot, the ratio of radial velocity dispersion to mean squared tangential velocity is constant. 



It has been conjectured (Ostriker & Peebles 1973) that the global stability of disks is related to the ratio of the 
total rotational energy to the total potential energy. As part of our aim is to examine this claim, let us derive the 
global virial quantities for future reference. The total kinetic energy K{R) and potential energy W{R) within any 
radius R can be computed as 



2 + 7 



K(R) = ttY^ovIR?,— 

^ °(l+7 + 2/3)(l-2^) \R^ 



R 



1-2/3 



W{R) 



2TTj:ovf,R'Q /B_ 
1-2/3 \Ro 



For the power-law disks, the virial theorem takes the form 

2 + 7 



2K{R) 



1 + 7 + 2/3 



W{R) = 0. 
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(35) 



(36) 



Note that the power-law disks do not, in general, satisfy the standard virial theorem 2K + W — Q. This is because 
it is not possible to "enclose" the system with a sufficiently large container. No matter how large the container, if 
the disk is warm, some stars will always cross its surface. When the disk is perfectly cold, the stars have no radial 
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motion and thus do not cross the surface. In this case, as seen from eq. (p6[), the standard virial theorem does hold. 
(Similar comments hold good for the isothemal sphere, for example). The rotational kinetic energy is 
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This concludes our summary of the equilibrium distribution functions for the power-law disks. In passing, let 
us emphasise that they are just the simplest choice of distribution functions, not the only ones. To illustrate this. 
Appendix ^ gives an example of a very different set of distribution functions. These perhaps fall into the class of 
remarkable curiosities, since the disks are built from orbits of one shape only. 



2.5 The Cut-Out Distribution Functions 



So far we have considered only the self-consistent case. We also plan to examine disks where parts of the central 
density are carved out. This is very much in the spirit of Zang's (1976) pioneering investigations. The cut-out mass 
is still present, in the sense that it contributes to the forces experienced by the remaining stars, but it is not free 
to participate in the perturbation. The disk is thus divided into "active" and "inactive" components. Although 
motivated partly by mathematical convenience, this is also a physically reasonable step to take. Stars in galactic 
disks are subject not merely to the disk's gravity field, but also to forces from the halo and bulge. A self-consistent 
distribution function, such as /s, is appropriate only when the disk's self-gravity overwhelms the gravitational 
potential of the other components. The immobile central mass can be interpreted physically as the hot bulge at 
the centre of disk galaxies. Another possibility - suggested to us by Tremaine (1997, private communication) - is 
to interpret the rigid density as caused by stars on highly elongated radial orbits. They pass through the centre 
of the disk, but they spend most of their time sufficiently far away from the disk so that they do not respond to 
the changing potential. There is one further motivation to consider cut-out disks as well as self-consistent ones. 
We wish to compare our results with those from iV-body studies - for example, those of Sellwood, Earn and 



collaborators (Sellwood & Athanassoula 1986; Earn 1993; Earn & Sellwood 1995). Obviously, numerical simulations 



cannot cope with infinite forces so in these studies the singularity at the origin is softened. Although some simulations 
(e.g., using Bessel function expansions) can track the behaviour of stars to infinite distances, disks with infinite mass 
are still problematic. For this reason, the disk is truncated at some finite radius. For comparison with these studies, 
an outer cut-out is needed, as well as an inner one. 



To ensure that we do not run afoul of Jeans' (1919) theorem, the carving-out is performed by multiplying our 
self-consistent distribution function fs (26) by a cut-out function H{Lz): 



fmtout{E,Lz) — H{Lz)fs{E, Lz). 



(40) 



We carry out our analysis for three cut-out functions H{Lz): the self-consistent (scale-free) disk, the inner cut-out 
disk and the doubly cut-out disk. These are given respectively by 
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(41) 



In Section 3, we will frequently use an equivalent form of the cut-out function expressed in dimensionless variables 
defined by H{Lz) = H{L). In all these formulae. 



Nb = 
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(42) 



where N and M are the inner and outer cut-out indices, respectively and must be positive integers. The choice of 
the inner cut-out reduces to Zang's (|l976| , eq. Z2.57) when (3-^0. Our generalisation seems to come "out of thin 
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air", but we shall see in Section 3.5 and Appendix O that our choice enables the analytic, rather than numerical. 



evaluation of a contour integral. Earn (1993) uses a doubly cut-out function of the same form as (41) to carry out 
his numerical simulations (although, not having to perform any contour integrals. Earn is free to choose non-integral 
N and M). 

What effect does this modification of the distribution function have on the surface density? When the disk is 
cold, the cut-out function depends only on radius and the active surface density may be calculated exactly. For a 
doubly cut-out disk, it is 
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(43) 



~ . . ~ "2/(2 — 13) 

where the truncation radius Rc is given by Rc — Lc ■ The proportion of the equilibrium disk which remains 
active rises from zero at the centre of the disk to unity at larger radii. At i? = i?o, the inner cut-out removes exactly 
half the equilbrium density. The value of N controls the steepness of the rise: the cut-out is gentler for smaller values 
of N. Conversely, the outer cut-out, parametrised by M, removes matter from the outer regions of the disk. At 
R = RcRq, it removes half the equilibrium density, with M controlling the sharpness of the cut-out. 

The active surface density of a hot disk can be calculated by numerical integration. In fact, heating the disk 
makes little difference to the active density. Figs. ^ and || show the active surface density of inner and doubly cut-out 
hot disks. In each case, the active surface density is expressed as a fraction of the equilibrium surface density. The 
form of these curves is close to that of H{Lz). Thus the value of H{Lz) is a good approximation to the proportion 
of density which is active at i? = RqLz- Our motivation for introducing the outer cut-off is to enable our results to 
be directly compared against A^-body work. In practice, the outer cut-out does not usually have a significant effect 
on the stability properties of the disk (unless it is so sharp as to provoke edge modes). 



3 THE FREDHOLM INTEGRAL EQUATION 

This section derives the Fredholm integral equation for the linear normal modes of the power-law disks. The first 
step is to decompose the imposed density into logarithmic spirals, as in Section 3.1. Linearising the collisionless 
Boltzmann equation enables us to calculate the density response. The condition for a self-consistent normal mode 
is that the imposed density equals the response density. This yields the Fredholm integral equation written down 
in Section 3.2. The kernel of the integral equation is the transfer function, describing how an excitation at any 
wavenumber feeds into response at all other wavenumbers. The general formulae for the transfer function is derived 
in Section 3.3. It depends on the properties of the equilibrium disk through the angular momentum function, as 
discussed in Section 3.4. 



3.1 The Logarithmic Spirals 

The logarithmic spirals were made famous by Snow ( |1952| ) and Kalnajs ( 1965| ; 1971). They have surface density 



/ f> \ 3/2-ia / P \ 3/2~ia 

Er™p = Epe^^e™^'-"^*) i^j = Epe*(™'^— *) , (44) 

where Sp is a constant amplitude. This m-lobed pattern rotates with a constant pattern speed ilp and grows or 
decays with a growth rate s. The logarithmic wavenumber (hereafter abbreviated to just wavenumber) a controls 
how tightly the spiral is wrapped. The pattern is more tightly wound for larger values of a. It is often convenient to 
collect the growth rate and pattern speed into the complex frequency to — mQp + is. By adding logarithmic spirals 
with different a, we can build up general patterns with azimuthal wavenumber m: 

/ + 00 
-oo 

Any imposed density perturbation Sjmp is expanded in a Fourier series of azimuthal harmonics characterised by 
m. In the linear regime, the response of the disk to the change in potential has the same order m of rotational 
symmetry as the imposed density perturbation. This means that our investigations can be confined to a single value 
of TO at a time. The expansion in logarithmic spirals is equivalent to taking the Fourier transform in the variable 
X = \ii{R/Rq). This is evident on defining 

= Spe^(™'^-"'*)e-3s/2 (46) 
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Figure 2. The active surface density of a liot disk with an inner cut-out function, shown as a fraction of the equilibrium density. The 
cut-out becomes more abrupt with increasing A^. [The solid lines are the density for A' = 1, 2, 3, 4. For fi = +0.5, = 0.199; for 
/3 = —0.5, CTu = 0.740. These are the temperatures at which the self-consistent disk is locally stable to axisymmetric disturbances.] 
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Figure 3. The active surface density of hot disks with a doubly cut-out function, shown as a fraction of the equilibrium density. The 
arrows indicate the direction of increasing M, or increasing severity of the truncation. The dashed line just visible at the top of the plots 
is the corresponding inner cut-out disk. [In each case, N = 2, Rc = 50. The solid lines are the density for M = 1, 2, 3, 4. The dotted line 
marks the position of Rc. For /3 = -1-0.25, 5-„ = 0.283; for /3 = -0.25,ct„ = 0.509.] 



SO that we obtain the conventional Fourier transform pair: 
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Here, Ajmp is the Fourier transform of the imposed density perturbation. For the transform to exist, Simp d) R^^^ 
must tend to zero as i? ^ and R ^ oo. This means that S 
R~^/'^ and must fall off at large radii more quickly than R~^/^ 
corresponding to a single logarithmic spiral component as 
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where K(a, m) is the Kalnajs gravity function 
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This function is real and positive for real a. It has the symmetry K(a,m) ~ K{~'a,m). 

From the imposed potent ial ip^^, we obtain the ch ange in the distribution function /j^. The linearised colli- 
sionless Boltzmann equation (Binney & Tremaine 1987, chap. 5) relates the change in the distribution function to 
the changes in energy and angular momentum experienced by individual stars as a result of the density perturbation: 



dt 



dtpu 



dR 



R 09 



9/ 
dE 



09 dL, 



df dE 
"dE~dt 



df dL, 
dL, dt 



(50) 
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In physical terms, eq. ( pO[ ) simply states that the rate at which stars move to perturbed orbits is equal to the rate 
at which stars leave equilibrium orbits. Integrating ( pO| ) over the entire time of the perturbation, we obtain 

fS^,{t)=-^AE-^^AL., (51) 

where AE and AL^ are the changes in the star's energy and angular momentum due to the perturbation. For a 
single logarithmic spiral component, these simplify to 

Physically, AE is the difference between the potential here and now, and an averaged potential sampled by the orbit 
over its history. Note that this derivation assumes that the perturbation vanished in the distant past. 



3.2 The Integral Equation 

To find the change in surface density S"o™ caused by a single logarithmic spiral component an integration of 

/j^p must be performed over all velocities u and v 

^?e7 ^ JJ fSZdudv. (53) 

To find the total change in density caused by the whole disturbance Eimp = J daAi^p (a) 'Sf^, we integrate over 
all the logarithmic spiral components 

daAi^p(a)E,"4". (54) 

-oo 

It is possible to equate Srcs and Ejmp and seek to solve the resulting equation for the self-consistent solutions. A 
more profitable approach is to equate the density transforms. We define the response density transform Arcs (o^) 
analogously to (47) so that 

daA.es(a)e^"^ A,,{a) = — di^e^-. (55) 

Substituting for Sics from (p4|), the response density transform becomes 

1 p+oo —iax r + oo 

A,,,ia)^— di—^ da'Amp(«')S?e's'". (56) 



-'m J —oo 



Exchanging the order of integration, we obtain (cf. Z3.42): 



+ 00 /' + 00 



1 ,^ E". 



Aos(a)= / da'^imp (a') (a, a') , where Sm{a,a') ^ — I dx -^e (57) 

Self-consistency requires the response density to equal the imposed density, or equivalently Ajmp (a) = Arcs (a)- 
Equation (57) therefore becomes an integral equation for the density transform of the normal modes of oscillation. 



Of course, similar integral equations have been derived by others before. For example, Kalnajs (|1971) derived 



a completely general integral equation in action-angle coordinates, whereas Palmer & Papaloizou ( |1990D res tricte d 



their attention to disks built from epicyclic orbits. This form of the integral equation was first derived by Zang ( 1976 ). 
The kernel of the integral equation Sm (a, ct') is called the transfer function and describes how much of the disks 
response to the disturbance with wavenumber a' occurs at the wavenumber a. It is to the transfer function that our 
attention now turns. 
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3.3 The Transfer Function 

This section contains an elaborate calculation of the transfer function of the power-law disks. To aid the reader, a 
reference table of some of the quantities are collected in Appendix The first step is to evaluate the change in 
distribution function given in eq. (^). From eqs. (^) and (52), the changes in energy and angular momentum 
caused by a single logarithmic spiral component are 

AE = i^Zlit) + tco f Cm" {t') dt\ AL, = im f ViTp (^0 dt' . (58) 

J —OO J —QO 

The integration is simplified by shifting to the frame which rotates at the star's average velocity ft. In this frame, the 
star's orbit closes and the integrands in eq. (58) become periodic with the radial period T. Following Zang (1976), 
let us define two dimensionless coordinates describing the periodic excursions of the stellar orbits as: 

X = In — = In ^, Y = 9 - nt = 6 - ni. (59) 

Rh 

Fig. ^ shows how X and Y vary over four radial periods for a somewhat eccentric orbit. The potential ■(/'imp ( |4^ ) 
experienced by the star at time t' when its position is {R', 6'), or equivalently X' = Ini?', Y' = 9' — fit' , is: 

i;^^^{t') = 27rGSpi^(a, m)RoR'^^^ exp \^i{mn - uj)t' + imY' + [ia ~ \)X'^ . (60) 

Now X and Y both have period T, so the term exp ^imY' + {ia — \)X'^ is similarly periodic as t' varies. It can 
therefore be expanded in a Fourier series (a method previously used by Kalnajs (1972) and Zang ( |1976D ): 

exp{imy' + (ia-i)l'}= ^ Qi™(a) exp i L (61) 

where the Fourier coefficient Qim{o^) is given by 

Qim{a)^^J^ cxp|^my' + (m-i)X'-^^'|d^'. (62) 



Qhn{a) - ^ / exp {zmf + {ia - ^)X' - tlx'} dx' ■ (63) 



Changing variables to the orbital phase x = I'^t, the Fourier coefficient becomes (Z3.27) 

2^ Jo 

Substituting this into our expression (|6^), we obtain for the perturbation potential sampled by the star 

4-00 

C^^(i')=27rGSpX(a,m)i?oi?^"^e'('"^'-"^*' Q;-(«)e*'"*'- (64) 

/— — OO 
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Using (58), the changes in the energy and the angular momentum of the star caused by a single logarithmic spiral is 



AE = 27rGJ:pK{a,m)RoRH 'e 



2 „i(mSl-Lj)t 



/ — — oo 



Ik + mfl 

Ik + mQ — oj ' 



(65) 



Ai, = 2^mGI]pi^(a,m)i?oi?t^"^e'(™"-")* V Jl^IIii^^f . 

^-^ Ik + mil — uj 

/— — oo 



(66) 



These correspond to (Z3.29) and (Z3.30), although Zang uses a slightly different form of the Fourier coefficient. 
Substituting these expressions into (pl|) gives the final expression for the change in the distribution function: 



/— — oo 



Ik + mil — uj 



dE dL, 



Thus far, our results are independent of the form of the equilibrium disk. Let us now specialise to the power-law 
disks by substituting for the derivatives of the distribution function (^) . Bearing in mind that the energy E has the 
opposite sign to /3, we can combine the derivatives obtained for positive and negative 



dE-~^ 



1 7 _ 7 
P P 2 



H{L,)L1\E\^/^+^l'^-^l^-^, ^ = CLl\E\^ 



1//3+7//3-7/2 



2 / H{L:,) dH 



We also replace the dimensional quantities with their dimensionless analogues (eqs. (12), (|22|), (|24)) 



7l, I ■ 



K — K L"^ . 



i?0 



dH _ 1 dH 
dLz VfjRo dLz 



Ri 







E^Ev}Li-\ 



~ Vf3 



(69) 
(70) 



Note that uj is defined slightly differently from k and fl. With the present definitions, all three dimensionless 
frequencies are independent of Lz (it is clear from eqs. (l2^), (|2^ ) and ( [l8|) that k and depend only on C7). With 
these substitutions, the change in the distribution function now becomes 



/-^(t) = 2^GSpX(a,m)Ci?J+^^ 
;=-oo + rnfl — ujLJ^ 



■|£.|l//3+7//3-7/2 



{Ik + mfl) 



1 7 _ 7 
/? 2 



7ml H(Lz) — mLz — — 

\E\ ^ J ^ 'dLz 



(71) 



This is the change in the distribution function brought about by a single logarithmic spiral component. 

Let us now proceed to find the transfer function Sm from (57). Substituting for the response surface density in 
terms of the response distribution function, we have: 



1 



/ + 00 ^ — iOLX 

dx— — 



ftn^pdudv. 



(72) 



This triple integral is transformed to one over the eccentric velocity U, the orbital phase x ^-nd the dimensionless 
angular momentum Lz- A careful calculation of the Jacobian (plead 1997) reveals: 



diS:,u,v) _ v$MU)U 

-6 Ljz 



d{x,U,Lz) 27r R 



(73) 



The transfer function then becomes 



Sm{a,a') ^ ^ dUdxdL^ ^ 



2^ ^Lr'ft^;\Lz,u,x)- 



(74) 



The angular momentum is integrated from zero to infinity, and the orbital phase from to 27r. The eccentric velocity 
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is integrated from zero to infinity for negative /3; for positive P the upper limit is (2//3 — 1)^/^. For brevity, these 
limits of the integration are not shown explicitly. 

We shall work in terms oi X = \nR = \n{R/RH) = x — ^t^Lz- Then R = e-^ , and on substituting for E,„ 
from (46) we obtain 



Sm{a,a') 



1 

27ri;p 2tt 



Note that t, X and describe where the star is in its orbit. They therefore depend on the orbital phase x, as well 
as on il which describes the shape of the orbit. Substituting for J^^^ from ([n]): 



'\ - f(i+7) GK(a',m)C 

27r 



Sm{a,a)^Rl w 



dL 



dU dx — — e 



± -im{e -at) -^iia+^)X ^ (a- a') In 



UMU)\E\-^' 



E 



i=~oo Ik + mf2 — u)Lz~'^ 



(Ik + mD,) 



1 7 _ 7 
f3 P 2 



— 7ml H(Lz) — mLy — — 

\E\ ^ i ^ 'dLz 



(76) 



Substituting E = {PU^ + /? - 2)/(2/3) and 9 - Qt = Y{x), 

Sr„{a,a') ^ Rl-^^vf^^^^-KGKia' ,m)C j dtjJo{U)U 



PU^+P-2 



+ 00 



^ Qim{a')— / dxe 



1 — — 00 

1 f dL, 



(Ik + mil) 



2/3 



2 + 27-^7 



^;72 + /? - 2 



(77) 



dH 

7TO H(Lz) - mLz^^ 
dL, 



From the definition of the Fourier coefficient (|6^), the integral over orbital phase x is just the complex conjugate of 
Qim{a)- Then, the transfer function becomes: 



S„^(a,a') ^ R2^'^vf/^^''^2TrGK{a',m)C / dUJo{U)U 



PU^+P-2 



2/3 



1 4- "J" _ 2. 



^ Qim{a)Q*i^^{a)Fhn{a~ a). 



I — — 00 



(78) 



Here, to make this expression a little more manageable, we have defined the integral over angular momentum to be 
the angular momentum function Fimir]), where r] — a — a': 



2^ 7o 



Ik + mfi — LoL, 



2 + 13 



[Ik + mfi) 



2 + 27 - /?7 



/3C/2 + /3 - 2 



1 ~ ~ ~ dH 

jm > H{Lz) - mLz^^ 
J dLz 



dLz 
L, 



(79) 



For convenience, the Fourier coefficients Qi„l and the angular momentum function _FJ,„ are shown as depending 
only upon the wavenumber a in ([78[), although they also depend upon the eccentric velocity U as well as on the 
disk parameters. The corresponding results for the transfer function and the angular momentum function for the 
Toomre-Zang disk are (Z3.40): 



iSm(a, a') 



7/2 



dUJoiU)Ue- 



Qun{a')Q*i^{a)Fir,,{ 



a ~ a 



(80) 



ZTT 



In 



Ik + m^l — ujLz 



{Ik + ■m9){l + 7) - 7m| H{Lz) - mLz^^ 
J dL, 



dLz 
Lz 



(81) 



Can we gain some rather more intuitive understanding of this complicated expression for 5m? The transfer function 
5,„(a, a' ) describes the contribution of the imposed logarithmic spiral component with wavenumber a' to the response 
component with wavenumber a. To see how this is calculated, we consider a star orbiting in the disk. The shape 
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Figure 5. Graph of the angular momentum function F22 plotted against 77 for two different growth rates in a /3 = 0.25 disk. In the left 
panel, s = 10~^; on the right, s = 0.3. In each case, the solid line is the real part, and the dotted line the imaginary. [Numerical details: 
P = 0.25, TV = 2, M = 6, Lc = 31.6; 7 = 13.7, Hp = 0.5, m = 2; I = 2, U = O.5.] 



of its orbit is characterised by its eccentric velocity U. The Fourier coefficient Qi^ describes the "match" between 
a particular logarithmic spiral component and this orbit. Specifically, the changes in the star's energy and angular 
momentum caused by the logarithmic spiral perturbation are expanded in harmonics of the orbit's radial period. 
The Fourier coefficient Qim gives the contribution to the ^th component of the perturbation experienced by the 
star. Thus, in the expression for the transfer function (|7^), the first Fourier coefficient Qim(a') describes how far 
the star is forced out of its unperturbed orbit by the imposed perturbation, and hence the star's tendency to stop 
contributing to the imposed logarithmic spiral. The second Fourier coefhcient Q/*„(q!) describes how well matched 
the perturbed star is to the response logarithmic spiral component with wavenumber a. Both these depend only on 
the shape of the star's orbit, not on its size; nor has any consideration yet been made of the "interaction" between 
the logarithmic spirals. This is accounted for by the angular momentum function Fi^- How easy it is for a star 
to move from the imposed to the response logarithmic spiral depends on the difference between the response and 
imposed wavenumbers, rj = a — a' , as well as the charactistics of the perturbation (its azimuthal symmetry, growth 
rate and pattern speed, given by m and to), and also on the size of the orbit, Lz- The integral in the angular 
momentum adds up similar-shaped orbits of all different sizes. Fim then describes how feasible it is for density to 
move from wavenumber a' to a in a perturbation with this m and uj. Owing to the scale- free nature of the disk, we 
have been able to deal simultaneously with all orbits of a given shape, irrespective of their size. The [/-dependent 
parts of the integrand in Sm measure how many stars there are for each shape of orbit. These are then added up to 
determine how much density ultimately moves from the logarithmic spiral component with wavenumber a' to that 
with wavenumber a. 



3.4 The Angular Momentum Function 

Our job is not quite yet complete! The angular momentum function eq. ( [79|) can be worked out analytically for the 
cut-out functions H{Lz) given in section |]^. This contour integration is by no means trivial and it seems wise to 
relegate the details of the calculation to Appendix In this section, we aim to understand the behaviour of the 
analytic expression for the angular momentum function in physical terms. Fig. |] compares the angular momentum 
function for two different growth rates s. The left-hand plot shows Fim{ri) for vanishing growth rate, while the 
right-hand one has s — 0.3. This figure exemplifies the trailing bias of these disks. Where Ik -\- mCl is positive and 
the eccentric velocity is low, the angular momentum function is highly asymmetric about = 0: the decay for 77 > 
is greatly attenuated. For negative values of lk-\-m(l or for high eccentric velocities, Fim decays rapidly as moves 
away from zero in either direction. The mathematical origin of this asymmetry is seen most simply in the expression 
for Fim in the case of a scale-free disk (C9). This depends on rj as 
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where fy is defined as 277/(2 + /?) (see Appendix y, (C4)). If the growth rate is zero, then on taking the principal 
value of the logarithm, we obtain 

F,™(7?)cxe rK + mn>0; (83) 

F,„(77) cxe """l^^l^ lR + mn<0. (84) 

1 — e^^^ 

Considering extreme values of 77 on either side of zero, it is apparent that for Ik + mil > the angular momentum 
function is be asymmetric about 77 = 0. The magnitude of Fim tends to a constant value for large positive 77, whereas 
it decays rapidly as exp(— 27r|77|) for negative 77. Conversely, for Ik + mfl < no such asymmetry is apparent; Fim 
decays as exp(— ttI?)!) in either direction. After the summation over radial harmonics and integration over eccentric 
velocity, the magnitude of Fim is likely, on average, to be greater for positive 77 (a > a') than for negative. This in 
turn means that, typically, the magnitude of the transfer function iSm(a, a') is greater for a > a', as demonstrated 
in fig. ^. Most of the response at a is due to imposed components with wavenumber less than a. An alternative 
way of viewing the situation is that imposed components mostly go to make up response components with larger 
wavenumber than their own. Thus the asymmetry in the transfer function means that our disk tends to make 
whatever pattern is imposed on it more trailing. 



3.5 The Special Case of Axisymmetry 

In the case of an axisymmetric perturbation, the integral equation admits certain symmetries. The first simplification 
is that the frequency is purely imaginary: lu = is. Secondly, the integral equation is equivalent to one with a Hermitian 
kernel, and hence must have purely real eigenvalues. It is straightforward to deduce the symmetries 

Fioiv) = Flu>i-V), Qro(") = 0«(-a), Q^io{a) ^ Q,o{a). (85) 

Let us now recast the integral equation (p5| ) into a form with a Hermitian kernel by defining a modified transfer 
function and density transform (Zang 1976) 



Ma, a') = J |^1^5o(a, a') B{a) = ^K{^)A{a). (86) 



The integral equation (|9^) now becomes 



n+oo 

XB{a) ^ I daB{a)Tf){a,a). (87) 

-00 



The value of this transformation is that To(a, a') is Hermitian, viz 

ro*(a', a) = To(a, a'). 
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The modified integral equation ( |87| ) is thus a homogeneous, hnear Fredholm integral equation with a Hermitian 
kernel. It must therefore have real eigenvalues and orthogonal eigenfunctions (Tricomi 1985). This occurs because 
the mathematical eigenvalues have no dependence on pattern speed. In general, the mathematical eigenvalue A is an 
analytic function of the complex frequency lj. Only for m = is the frequency purely imaginary and the problem 
one-dimensional. 

There is a further symmetry. We make the transformations a — > —a, a' —a', / —I, so as to obtain 
Tq(— a, — a'). Using the symmetry properties of the Fourier coefficients and the angular momentum function ( |85| ) 
along with that of the Kalnajs gravity factor, we obtain the result 



To{~a, ~a') = To{a, a'). 



(89) 



This has the consequence that the eigenvalues are degenerate. If B{a) is an eigenvector with eigenvalue A, then so 
is B*{—a). To see this, we make the transformations a 
take the complex conjugate. We obtain: 



-a, a' — > —a' in the modified integral equation (87), and 



XB*{-a) 



da'B*{-a')T^{-a,-a') = 



+ 00 



da'B*{-a')TQ{a,a'). 



(90) 



In other words, B*{—a) also satisfies the integral equation (p7|). In terms of our original density tr ansforms, A{a) 
and A* (—a) are a degenerate pair with the same eigenvalue A. This is just the anti-spiral theorem (Lyndcn-BcU & 



Ostri cer 1967; Kalnajs 1971) for axisymmetric perturbations. 



4 NUMERICAL METHODS 

This section discusses the numerical algorithms required to find the normal modes. Sections 4.1 and 4.2 consider 
the evaluation of the Fourier coefficients and the transfer function respectively, whereas Section 4.3 presents the 
discretisation of the integral equation. Again, let us emphasise that the numerical method is adopted - with a little 
streamlining - from Zang ( |1976D . 

4.1 The Fourier Coefficient 

The integrand in the definition of the Fourier coeffients ( |63|) is periodic with period 27r in x- We are free to define 
X = to correspond to pericentre. This means that at the time i = 0, the star has radial coordinate R — i?min- As 
X is even, and Y is odd, about x = 0, we can write 

Qiniia) ^ ^ cxp|(ia- i)A|cos(my-/x)c^X- (91) 



The equations of motion (16) are solved by fourth-order Runge-Kutta integration (Press et al. 198£, chap. 15) to 



obtain the stellar position as a function of time, and thus X and K as a function of x- The integration over x is 



carried out by the midpoint method (Press et al. 1989, chap. 4), using points in the midpoint integration and 
in the Runge-Kutta. The problem is to choose large enough to obtain excellent accuracy, while keeping it as 
small as possible in order to save time. Eccentric orbits need more work to obtain good accuracy, so, as suggested 



by Zang (|T97§), is made to depend exponentially on U : 

r^i/. = aaccexp(&acct/)- (92) 

Values of Cacc ~ 10 and 6acc — 1-5 usually worked well. With these values, the Fourier coefficient is obtained with 
around 6 s.f. accuracy for low eccentricity orbits (U = 0.5), but perhaps only 1 or 2 s.f. for high eccentricity orbits like 
{U = 1.5). The Fourier coefficients are generally smaller for higher values of the eccentric velocity. In the integrand 
for the transfer function, two Fourier coefficients are multiplied together at every eccentric velocity. The accuracy 
with which the transfer function is obtained is thus dominated by the accuracy of the Fourier coefficients at low 
eccentric velocities. This is useful, since the Fourier coefficients are costly to evaluate when the eccentric velocity 
is high. Fig. ^ shows the Fourier coefficient as a function of wavenumber, for two different eccentric velocities. It is 
shown only for positive wavenumber. Its behaviour for negative wavenumber is readily deduced from the symmetry 
property Qim{oi) = Qlmi—c^)- We see that U controls the amplitude and the frequency of oscillation of Qim- Fig. |^ 
compares Fourier coefficients at different radial harmonics. For high values of /, Qim remains close to zero until a 
is large; thereafter it oscillates with lower frequency. Logarithmic spiral components must be tightly-wound in order 
to excite responses at high radial harmonics. 
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Figure 8. Graph of the Fourier coefficient Q;,„ plotted against wavenumber a for four different radial harmonics I. In each case, the solid 
line is the real, and the dotted line the imaginary, part of Qim- [Numerical details: (3 = 0.25, m = 2, U = 0.5; aacc = 100, bacc = 2.5.] 



4.2 The Transfer Function 

The expression for the transfer function ( [78| ) involves a sum over radial harmonic / from — oo to +cxd. Fortunately, the 
magnitude of the terms in this sum decreases sharply with |Z|. Negative values of I decrease their contribution faster 
than positive /. A good approximation is given by summing from I = /min to I = Imax, where Zmin is negative and |/min| 
is less than Imax- Values of Imin = —20, Imax = 30 usually worked well. The transfer function also contains an integral 
over eccentric velocity U, which adds up orbits of all possible shapes. The e quival ent expression for the Toomre-Zang 
disk ( ^0|) contains a Gaussian factor exp{~ ^tj"^ / af^) . This prompted Zang (197£) to use Gauss-Laguerre quadrature 
to evaluate Sm ■ The fundamental formula is 







f{x)e-''dx^}^wj{x,), (93) 

1=1 

where f{x) is a smooth polynomial- like function. The weights Wi and abscissae Xi are well-known (see Abramowitz 



& Stegun 1989; Press et al. 1989, chap. 4). Zang ( p. 9 761) found it preferable to reduce the dispersion of the Gaussian 
in this expression, replacing it with (t„ = fa-d'u, where the fraction f„ is around 80%. This concentrates attention on 
the lower values of U, where the integrand is larger. For the power-law disks, the distribution of radial velocities at 
any spot is similar - but not exactly equal to - a Gaussian with dispersion (7„. This suggests that Gauss-Laguerre 
quadrature may still work well. Although our integrand is not directly of form (p3|), it is readily made so. The limits 
on our integral are zero and infinity for negative /3. For positive /3, the upper limit is U = (2//3 — 1)^/^. At least for 
the models of most interest in galactic astronomy (|/3| < 0.5), this distinction is in practice unimportant, since the 
integrand has fallen to negligibly small values well before this limit on the eccentric velocity. 

Let Jjj be the integrand in Sm when the integration is carried out over U, so Sm (a, a') cc JudU . We then 
have 

CT^ ( U'^\ 

X.=X,^exp , (94) 
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where Sm {ot, a' ) cc ^^Xye-'^d V and V ^ \U^/ al. The rehabihty of Gauss-La guerre quadrature was tested against 
the extended midpoint method (Press et al. 1989, eq. (4.1.19)). These comparisons (Read 1997) indicate that Gauss- 
Laguerre quadrature is a remarkably efhcient way of performing the integral. Excellent six significant figure accuracy 
is obtained with a handful of function evaluations. 



4.3 Computational Solution of the Integral Equation 

Given that the kernel of the integral equation can be evaluated to high accuracy, we now need to devise a method 
for its numerical solution. We ultimately seek self-consistent solutions, for which Aics{ol) = Aiinp(a). However, as 
Zang (1976) already noted, the first step is to consider the more general problem for which the response density is 
a complex multiple of the imposed density, i.e., Ares (a) = AAinip(a). This casts the integral equation (57) into the 
form 

p+ca 

XA{a)^ da'A{a')S„^{a,a'), (95) 



where we have dropped the subscripts on A. We refer to A as the mathematical eigenvalue. Of course, only the 
instance when A is unity carries the physical significance of a mode. The advantage of this mathematical artifice is 
that the integral equation ( p5|) is in the standard form of a homogeneous, linear, Fredholm equation of the second 
kind (see e.g. Courant & Hilbert 1953; Delves & Mohamed 1985). For a given value of A, such an equation normally 
admits only the trivial solution, A(a) = 0. The values of A for which non-trivial solutions exist are the eigenvalues 
of the equation. We seek an iterative scheme which drives the eigenvalue to unity, thus providing a self-consistent 
mode. There is a close analogy between linear algebraic equations and linear integral equations. The former define 
relations between vectors in a finite-dimensional vector space, the latter define relations between functions in an 
infinite-dimensional vector space (technically, a Banach space). This analogy can be made explicit by applying a 
quadrature rule to the integration in ( psf ) to obtain 

XA (aj) = ^ WiA (ofj) 5„i(aj, tti), (96) 



where Wj ar e some appropriate weights. This general approach is called the Nystrom method ( |Delves fc Mo- 
hamed 1985| ). It evidently reduces the solution of an integral equation to the solution of an algebraic eigenvalue 



problem. The latter is a classic and well-studied area of numerical analysis, for which many tried and tested tech- 
niques are available. 

Much of the skill in the numerical solution of integral equations comes from the choice of quadrature rules and 



weights. Zang (1976) devised an elegant method based on locally approximating the kernel and the response by 
Lagrangian interpolating polynomials. This is naturally adapted to the instance when the kernel varies on a much 
smaller scale than the solution. We follow this method, but we did not find Zang's splitting of the kernel of the 
integral equation into Hermitian and Volterra parts to be needed. First, to obtain smoother functions, the Kalnajs 
gravity factor is extracted by defining 

Smia,a') = K {a\m)Sm{a,a') , A{a) = (97) 

K (a, m) 

The eigenvalue equation now becomes 

/-I-CX3 
da'A{a')S„^ia,a') . (98) 
-oo 

Let us introduce a finite grid of points in wavenumber space, a^. If we need to know the value of A or Sm at a value 
of a intermediate between the gridpoints and a^+i, we interpolate over the eight gridpoints from ars to Ur+A- 
For 8 equally-spaced points Aa apart, Lagrange's classic formula for the interpolating polynomial P{a) through N 



points f{ak) (Press et al. 1989, chap. 3) becomes 

8 



P{a) = n (a - a.+.-4) ^3 - k)l iAa)^a-lr'^- kAa) ' ^^'^ 



Defining x — (a — ar)/ Aa, this becomes 



8 



P{a) + E^(-lf (3^,),{frfc)l(^_,) - E^i.N/(a.+.), (100) 
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where 



nil (x + 4-^) 



{3 + ky.{A-ky.{x~k)' 

Then the interpolated approximations to the response and the transfer function arc (cf. Zang 1976, app. D) 

4 4 



(101) 



i(a') « Lk[x']A{ar+k), 



(102) 



fe=-3 



k=-3 



where x' = (a' — ar)/Aa. The infinite range in the integral equation ( |95| ) is a kind of singularity. Truncation of the 
wavenumber range at large finite values is the simplest but most brutal way of handling the singularity. In practice, 
we found this to be surprisingly effective. Wavenumber space is approximated by a finite grid with n points along 
each side. We choose n, and the grid-point spacing Aa, so as to obtain a sufficiently accurate solution of the integral 
equation (57). The finite size of the grid means that we encounter problems when evaluating ^i,np and Sm for values 
of a near the edges of the grid. We are interpolating over eight points, so we need data from a-2 to a„+3. We deal 
with this problem by simply assuming that our function is zero outside the grid. This is acceptable provided the grid 
is large enough. Then the values of the kernel at the missing interpolation points are negligibly small. This procedure 
is justifie d empiricall y by the demonstration that the mathematical eigenvalue converges to a value independent of 
grid size ( Read 1997 ). 

Substituting the Lagrange-interpolated approximations for Ai,„p and Sm (102) into the modified integral equa- 
tion (pa), we obtain 



4 4 

XA{a)=K{a,m) ^ ^ 



da' Li 



a — ttr 



Lk 



a — ttr 



i(a^+i)5m(a,a;+fc). 



(103) 



Of course, the integration over wavenumber now runs from a = ai to «„, instead of from a = —oo to oo. We break 
this integration into n portions of Aa. 



XA (a) ^ K{a,m)J2Yl / 
Then, changing variables to x' = (a' — ar)/Aa, we obtain 



Aa 



Aa 



i(a;+,)5„,(a,a;+fc). (104) 



n 4 4 „i 

AA (a) = X (a, m) ^ ^ A{ar+t)Srn{a,ctr+k)Aa dx' Li[x']Lk[x'] 

r=li=~3fc=-3 •'^ 



Let US define the weighting coefficients Cik by (cf. Zang (197G), app. D) 

Cjfc = Aa / dxLi[x]Lk[x]. 
Jo 



(105) 



(106) 



Although there appear to be 64 Cik, only 20 of them are independent, because of the symmetry properties Cik = Cki, 
and Cik ^ C'(i_i)(i_fc) . The weighting coefficient s are easily evalu ated using the midpoint method given in Press et 
al. ( 1989 , chap. 4). They are tabulated in Read ( 1997 ). Equation (102) can now be written as: 



XA (a,) = K (a„ m) ^^^^"^ 

{aj , ar+k)A{ar+i) 



(107) 



r=l j=-3 A;=-3 



Terms in this multiple sum for which r + i is less than 1 or greater than n contribute nothing. We can thus write the 
above equation in terms of a single sum from ai to a„, and collect the weighting coefficients Cik, smoothed transfer 
function Sm and density profile A together into a single quantity S. 



^Maj)=Y,SjsA{as)- 



(108) 



This is of course just a matrix equation for the mathematical eigenvalue A. If there is an eigenvalue of the matrix 
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S equal to unity, then the corresponding eigenvector A gives us a self-consistent mode, i.e. a self-sustaining density 
perturbation. Sjs cannot be simply expressed in terms of Cik, Sm and A. But we see that it represents the total 
coefficient of A{as) in eq. ( |l07l ). 

Having obtained the matrix Sjs^ it is reasonably straightforward to find its eigenvalues using the eigensystem 
package {EISPACK) developed by Smith et al. (1972). This returns the eigenvalues of a matrix and, optionally, the 
corresponding eigenvectors. The power method can also be used to provide an independent check on the largest 
eigenvalue. This relies on repeated application of the matrix S to a test vector t (Ralston & Rabinowitz 1978, chap. 
8, sec. 7). Probably the most important factor in achieving excellent accuracy in the eigenvalues is the size and 
fineness of the grid. The range of the grid, (n — 1)Aq;, typically needs to be around 50 in order for the mathematical 
eigenvalue to be accurate to 6 s.f. A grid-spacing of Aa = 0.2 is usually sufficient. However, suitable values necessarily 
depend to some extent on the form of the mode. For instance, for high azimuthal harmonics, the modes tend to be 
tightly-wrapped, requiring a more extensive grid. 

Having established that the mathematical eigenvalue can be found to good accuracy, let us consider how to locate 
the unit eigenvalues which correspond to self-consistent modes. We are investigating modes of a given azimuthal 
symmetry in a disk with given density profile. This means that m, N , M and Lc are held fixed. It leaves us free to 
adjust 7, s and Sip. A disk at a given temperature, if it admits a mode at all, generally does so only for a particular 
growth rate and pattern speed. We therefore hold 7 fixed, and adjust s and Sip until a mode is found. This is in fact 
a search for the root of the equation A = 1 in only one dimension, since the mathematical eigenvalue A is an analytic 
function of the complex frequency w = mfip -I- is. The Newton-Raphson method is extremely effective at locating 
such roots. When the disk is sufficiently hot, it is generally completely stable. As the disk is cooled, instabilities set in 
through the marginal modes, for which the growth rate s is zero. We are therefore often interested in the marginally 
stable modes. To find these, s is set to some vanishingly small value, and a two-dimensional search is performed for 
the critical pattern speed fip and temperature 7 using the Newton-Raphson method in two dimensions. 



5 CONCLUSIONS 

This paper has set up the machinery for an investigation into the large-scale, global, modes of the power-law disks. 
The Fredholm integral equation for the normal modes has been derived. This is done by linearising the coUisionless 
Boltzmann equation to find the response density corresponding to any imposed density and potential. The normal 
modes are given by equating the imposed density to the response density. This scheme is implemented in Fourier 
space, so that the Fredholm integral equation relates the transform of the imposed density to the transform of the 
response density. Numerical strategies are given to discretise the integral equation to yield a matrix equation. This 
requires some care as the kernel of the integral equation is either singular (as in the case of the self-consistent disk) or 
almost so (the cut-out disks) . The problem of locating the normal modes is thus reduced to an algebraic eigenvalue 
problem, for which a wealth of standard techniques are available. The crucial simplification underlying the analysis 
is that the force law describing the gravity field of the galactic disk is scale-free. It is this that ensures that orbits 
passing through any one spot in the disk are simply scaled copies of the orbits passing through any other spot. When 
the orbits are forced by a scale-free logarithmic spiral perturbation, the response of all orbits of the same eccentricity 
can be computed at the same time. It only remains to add up the contributions at every eccentricity. For a disk with 
an arbitrary force law, such simplification is not possible. 

Let us again emphasise that nearly all the computational techniques required to study the stability of the power- 



law disks were developed by Zang (1976) in his pioneering study of the disk with a completely fiat rotation curve. 
Our contribution consists only of deriving the extensions to study the complete family of scale-free power-law disks. 
This paper has discussed just the algorithms - the aim has been to present all the mathematical and computational 
details for reference here. The following paper in this issue of Monthly Notices implements the machinery discussed 
here to provide a complete description of the spiral modes of the self-consistent and cut-out power-law disks. A fast 
computer code is developed which gives the growth rates and pattern speeds of the normal modes for any power-law 
disk in a matter of minutes when run on a modern workstation. 
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APPENDIX A: SINGLE-ECCENTRICITY DISTRIBUTION FUNCTIONS 

In the main body of the paper, we have used stellar distribution functions for the power-law disks that arc built 
from powers of energy and angular momentum. Many other equilibrium distribution functions are possible. In this 
Appendix, we consider disks in which all the stars have the same shape of orbit, characterised by an eccentric velocity 
U = Ur- Mathematically, we look for distribution functions of the form 

U = F{R^)5{U -Ur). (Al) 

The surface density is 



^eq — So 



R 

We transform this to an integral over U and Rh, using the Jacobian 

dudv = dUdR, fl - f ) ^ i^)"' , . (A3) 

and substitute our assumed distribution function 

^ /SoX'*" (, /3^ 1-g ll(R«\"'- UF{Ru)S(U-U,:jdUdR„ 
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We integrate over U , and then transform the integral over _Rh to one over R = R/Rh- 

•^fln.„ R-^Ju^ + l-R-^ + ^ARP -I) 



J V 2 

where Ur is the value of Ur in diniensionless units. Looking at the powers of R , we can guess at a solution of this 



integral equation, namely F(R/R) — k{R/ Ry^^^'^ . Substituting this into (A5), and remembering the definition of 
the auxiliary integral Jn{U) (|l8|), we can solve for k and obtain 



1+/3/2 ^ pl+/3 

(A6) 



The distribution function is then 



MU,Ru) = - ^^^^^2 7^ (A7) 



This expression is valid for /3 = 0, in which case the appropriate form for the auxiliary integral J7i must be 
used. 

All the stars in this disk have orbits of the same shape, but different size. Each star sweeps out an annulus as 
it orbits. The relative width of this annulus depends on the eccentric velocity Ur, but the overall size of the annulus 
depends on i?H- It is easy to imagine that if we add up annuli with a variety of i?H in suitable proportions, we could 
recover the surface density of the equilibrium disk; and this is what the distribution function ( [A7| ) does. It seems 
intuitively right that the number of annuli needed falls off with i?H at the same rate as the surface density falls off 
with R. Similarly, we can see why the total number of annuli depends inversely on eccentric velocity Ur'. for high 
Ur, the annuli are very wide, and few are needed; for low Ur, the converse is true. 



APPENDIX B: REFERENCE TABLE OF DIMENSIONLESS QUANTITIES 



Quantity 



For the Toomre-Zang disk 



For the general power-law disk 



-Rh Home radius 

i Time 

u Radial velocity 

V Tangential velocity 

E Energy 

Lz Angular momentum 

C/ Eccentric velocity 



Rh 

t = 

u = 
V — 

E = 
Lz = 
U = 
U^ 



Rh/Rq 



Ru 
dR _ 
dt 

Rde 
dt 

E . 



R- 



E : 







In 



Rovo 
_ U_ 

= -S^ - 1 



= R 



R- 



2\uR 





Auxiliary integral 




2dR 


fl"((72 + l-_R-2 








Radial frequency 


k = = 


2ir 
Jo 


n 


Angular frequency 


n = = 

Vo 


Jl 
Jo 




Epicyclic frequency 


ko — y/2 




no 


Circular frequency 






X 


Orbital phase 


X — Kt — kt 




Y 


Angular deviation 






X 


Logarithmic radius 






X 


Scaled logarithmic radius 




) = InR 




Radial velocity dispersion 


5-2 = ^ 

»' 1+7 




1 


Anisotropy parameter 







Rh 
t = 
u — 

V — 

E = 

Lz 
U 

k = 



— Ru/ Ro 

dR. u 5/3/2 

IF = 

Rde 



Vfi ' 



''H 



HdU V 



H 



R- 



R-H ' LJ 



Rovp 
^U^ -I 

flmax 



R-2 2 

R - p 



[R"0 

2dR 



1 



1 



R"'{U'^ + l-R-'^ + f^{R-l3-l))^ 
Rh K)/5/2^ _ 27T 
Jo 



''H 



Jo 



Clo 

X = 
Y = 

X = 
X = 



= 

= 1 

' Kt ~ kt 



In 



1+7+2/3 



Ini? 



1 - 2/3 
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APPENDIX C: THE ANGULAR MOMENTUM FUNCTION 

In this Appendix, we derive the forms of the angular momentum function by explicitly performing the contour 
integration over angular momentum L^- Let us recall the angular momentum function is 



1 



Ik + mfl — LoLz 



(Ik + mQ) 



2 + 27-/37 



l3U^+P-2 



7771 ( n — mLz — — 
dLz 



dLz 
I. 



The integration is carried out over the variable h defined by 



H{Lz)^H{h). 



(CI) 



It is performed for three different cut-out factors, namely the self-consistent disk, the inner cut-out disk and the 
doubly cut-out disk (all defined in eq. (41)). We proceed by splitting up the integral into two parts: 



Firniri) <{lfi + mCl) 



2 -I- 27 - /37 



7m > Ji — mJ2, 



(C2) 



where 



12-/3 e^i-n2h'' - 1 

J,^f3^r,,L,) = —--L n{h)dh, J,^f3^r,,L,) = — 

2tt 2 + fj Ik + mfi - use'^ 2tt 



dn 



Ik + niQ - ue^ dh 

We need evaluate these only for the case (3 = 0, since the general integrals are related to the /3 = integrals as 

2 



dh. (C3) 



2-/3 / ~^±M- 



J2(/3,ry,L,) = ( 0, r}, if^ 



with f] — rj X 



2 + /3 



(C4) 



Lc has been included as an argument, although clearly this is only relevant to the doubly cut-out disk. We suppress 
the third argument unless needed. 



CI The self-consistent disk 

Here 



1 r°° 

2t^ J-oo Ik 



-dh. 



J2(0,?;) = 0. 



(C5) 



We first consider the case ^ 0. Ji may be evaluated by contour integration. The integrand has poles when 
hn = ln[(ZK -|- m^)/Cj\ + 2n«7r, where n is an integer. The poles occur at intervals along a line parallel to the 
imaginary axis. Note that w is not in general real (w = mJlp -|- is), so that the poles are displaced from the lines 
h = 2niTT. We evaluate the integral 



1 

2^ 



^Xh^-irih 



Ik + mQ — uje^ 



dh. 



(C6) 



around a rectangular contour as shown in the left-hand plot of fig. CI, with long sides at ft, = 0, 7i = 2z7r, and short 



edges at ±L, enclosing the pole at ft = lii[{lK + mil) /uj]. Here the logarithm denotes the principal value, i.e. the 
imaginary part lies between and 2tt. Then 



1 

27r 7_L Ik + mn - uje^ 2ti 
1 f^'' e-^^e2"'e*^fe*''^dft 



_L Ik + mfl — uje'^ 



1 

2^ 



Ik + mil — uje^ye^ 



(C7) 



X In '- 



-ri In - 



Ik -f mfi — Cbe^ye"'^ Ik + mfi 
Taking L ^ 00, and then A ^ 0, we obtain the result 

Jiio,v) = - 



-in In 



{Ik + mVl) (1 - e^'T'J) 



(C8) 
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Figure CI. Contours and poles for the self-consistent disk. The left-hand plot shows the contour used for r) ^ 0, and the right-hand 
that for r] = 0. 



Generalising this result to all (3 as described above, we obtain 



2-/3 
2 + (3 



2 + 2j- (3^ 



(3U^+(3-2 



7771 



-ifjln - 



(C9) 



This expression is singular when 7/ = 0. As we now demonstrate, Fim{ri) contains a delta- function at r/ = 0. For the 
case rip = 0, s = 0, we can integrate the expression for Fim{tf) directly: 



^i(0,?7) = - 



1 



1 



Ik -h mVl 27r 



1 



Ik + mfl 



and so for the general power-law disk 

FiM = 



2-/3 

2 + /3 



2 + 27-/37 



/3f/2 + /3 - 2 



7777 



<5(7)). 



To obtain the expression for general ilp and s, we write 



Jl(0,77) = Ji(^\0,7/)-f ,5(7?), 



(CIO) 



(Cll) 



(C12) 



where J^^^ (0, 77) is equal to our previous expression del) for Ji(0, 7;). We integrate ( |cT^) over 77, thus smoothing out 
the delta-function. 



driJi{0,ri) 



77=-e 



d77 4'^(0,77) + A' 



(2) 



))=-£ 



4-c 



dr]5{ri), 



(C13) 



77=-e 



We now let e ^ 0. We see from the expression for j\'^\ rf) al ready obtained (|C^ ) that Ji^\ri) is odd in the limit of 
vanishingly small 77. It therefore contributes nothing in ( |C13| ). Using the expression (C3) for Ji(0,77), we obtain 



= lim- 



drj 



dh- 



-irjh 



Jh=-oo Ik + mVl - uje^ 
Swapping the order of integration and performing the integration over 77, we obtain 



= lim — 



1 



dh e'^'' - e""'* 



27r7 .lh=-oo h lk + - ue^ 



(C14) 



(C15) 



We evaluate this by carrying out two separate contour integrals. For the integral with e""'', we close the contour in 
the lower half-plane, as shown in the left-hand plot of fig. CI. As usual, we take the radius of the inner semi-circle 
to zero, and the that of the outer semi-circle to infinity. The integrand then vanishes on the outer semi-circle, and 
we obtain 



1 



dh 



1 



1 



1 



27r7 Jh=-oo h lk + mfl - uie^ Ik + mfl In IBd^ + 2m7r 2 Ik + mfl - Co 



(C16) 
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The integral with e+"^'' is evaluated over an analogous contour in the upper half-plane. We obtain 



dh 



1 °o „ie In i 



27r« h=-oo h lk + mQ - Cbe'^ Ik + mQ In Ih+mEl + 2m7r 2 + mD. - i 



Subtracting these two results and taking the limit e ^ 0, we obtain 



dh e'^'' - e-*"'' 



1 



1 

27ri Jh=-oo lk + mQ — uje^ In + mil - Co Ik + mQ 2mT — i In 



E 



The sum can be explicitly evaluated 



E 



n— — oo 



2n7r — iz 



E 



E 



2n7r — iz ^-^ 2n7r — iz z 

n— — 1 n— 1 



2*^E 



1 



4^ 



(n + iz/2'K){n — izjlix) 



1 iz 
- t: cot — 

2 2 



(2n7r)^ + 

i + 1 
2e-^ - 1 



where we have used the standard result ( [Prudnikov et al. 1986| , eq. (5.1.6.4)) 



E 



n— — oo 



(n + a)(n + b) b — a 



(cot 7ra — cot nb) 



Thus we obtain 



(2) 



1 



2(;k + mn) 

The angular momentum function for the self-consistent disk is then 



2-/3 



2 + f3 



2 + 27-^7 



/3C/2 + /3 - 2 



i^t;;;^; 2 — 



^—if} In - 



(C17) 



(C18) 



(C19) 



(C20) 



(C21) 



(C22) 

The corresponding result for the Toomre-Zang disk is (correcting a typographical error in Zang (1976), eq. (Z3.43)) 

(C23) 



■j — irj In - 



C2 The inner cut-out disk 

Here, we must evaluate 



1 r°° 

27r J_oo Ik 



r,Nh 



+ mn- Cjef^ e^f^ + 1 



dh. 



1 

27r J_oo In 



Ne 



Nh 



:dh. (C24) 



Again the integrands have poles when h = ln[(^K + ni^)/Cj\ + 2niTT. However, they now also have poles along the 
imaginary axis, at hj — (2j — l)i7r/iV, where j is an integer. We integrate around the same rectangular contour as 
for the self-consistent disk (fig. CI). As well as the pole at ln[(ZK -I- mfl)/uj], we now also enclose N poles lying on 
the imaginary axis. Adding up the residues from all these poles, we obtain the result: 



Ji(0,ry) 



(Ik + mQ) 



N-l -ir/ln ■ 



N 



e « 



(Ik + mn)^ + Cj^ 



E- 

,=1 ^K + mfi — we" 



(C25) 



^2(0,7;) = - 



{lk + mh)^-^NCj^e~'-'^^'^'- 



1 - e^^^'' ■ 



\{lk + m^Y 



N 



e « 



Ik + mVt — uje' 



IT] 



~. ~ ~ 2 j — 1 ■ 

Ik + mf2 — uje^^'' 



(C26) 
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For /3 ^ 0, we obtain the angular momentum function 



1 - e'^^^ 
2-/3 



N 



2 + (3 



{Ik + niQ) 



{Ik + mn)'^ + Cj^ 
2 + 27-/37 



/3C/2 + /3 - 2 



7TO 



N ^ IT. 



e « 



Trr; 



2j-l 



Ik + mfi — (2;e ' 



, 2-/3 



2 + /3 



(Zk + 777i7) 



2 + 27-/37 



/3C/2 + /3 - 2 



77)7, 



(C27) 



We can use I'Hopital's rule to obtain the result for rj = Q. This requires that both the numerator and denominator 



of (C27) are zero in the limit ?/ ^ 0. The denominator (1 — e?"^"^) is obviously zero in this limit, and it can be shown 
that the numerator is also (e.g. it is readily apparent for N = 1). Application of I'Hopital's rule then yields 



-F/™(ry = 0) 



i J i{lk + mn)'^~^\n 



27rl {Ik + mn)^ + 



N 



2-/3 
2 + /3 



{Ik + 777f2) 
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C3 The doubly cut-out disk 

We briefly summarise the results for the doubly cut-out disks. Here, 
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The integrand now has additional poles at hk — {2k — l)77r/M + (2 + f3)/{2 — f3) x InLc, where k is an integer. We 
thus obtain a second sum, and Fim becomes 
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In the limit ry ^ 0: 
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